/**
 * @brief   : 
 * @version : 0.1
 * @author  : wanghuan
 * @date    : 2024.10.26
 * @copyright Copyright (c) 2024 by 神工仿真, All Rights Reserved.
 */

#include "ThermalSolver.h"
#include <algorithm>
#include <iostream>

#include "Utility/Algebra/Factory.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"

#include "DBManager/IService/Result/Nodal/ITEMPERATUREService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/DataAdapter.h"
// #include "SGFem/ElementCalculator/Factory/ElementCalculatorFactory.h"
#include "SGFem/ElementCalculator/Register/ElementCalculatorRegister.h"
#include "SGFem/Task/Kernel/MModel.h"



using SG::DataStructure::Common::Status;
using SG::DataStructure::FEM::Project;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::FEM::Subcase;
// using SG::DataStructure::FEM::SPCADDData;
using SG::DataStructure::FEM::LoadCombination;
// using SG::DataStructure::FEM::ConstraintType;

using namespace SG::DBManager;

using SG::FEMSolver::EigenSparse;
using SG::FEMSolver::EigenVector;
using SG::FEMSolver::AssemblyEigen;
// using SG::FEMSolver::Element::ElementCalculatorFactory;
using SG::FEMSolver::Element::heatElementCalculatorRegister;
using SG::DataStructure::FEM::SOLType;



/// @brief static member
std::shared_ptr<SG::Algebra::LinearSolver> ThermalSolver::s_solver;

// SG::DataStructure::Common::Status ThermalSolver::CalcMatrixEquation (const DBServiceFactorySP dbServiceFactory,  const SG::DataStructure::FEM::Subcase& subcase, const std::set<SG::DataStructure::Common::Id>& loadCases)
// {
//     // 对应工况的分析类型必须是 `传热`
//     // 即 AnalysisTypes::ANALYSISTYPE_Heat == subcase.m_sol


//     // std::cout << "===== compute stiffness \n";
//     // const auto elementTypes = dbServiceFactory->get<IElementSummaryService> ()->GetAllElementTypes ();
//     // for (const auto& type : elementTypes)
//     // {
//     //     //此处跳过边界单元,后面根据边界条件进行计算
//     //     if ()
//     //     {
//     //         continue;
//     //     }
//     //     if ()
//     //     {
//     //         continue;
//     //     }

//     //     auto eleCalculator = ElementCalculatorFactory::Create (type);
//     //     if (NULL == eleCalculator)
//     //     {
//     //         std::cout << "can 't find element calc "<< std::endl;
//     //     }
        
//     //     eleCalculator -> Initialize (type);
//     //     eleCalculator -> ComputeElementMatries (dbServiceFactory);  
//     // }
//     // std::cout << "===== compute stiffness contribution finished\n";
//     // std::cout << "===== compute loads contribution finished\n";

//     // // 计算各个载荷集各自集中热源对载荷端的贡献,热源，固定热流
//     // for (auto& loadCase : loadCases)
//     // {
//     //     AnalysisProcess::ComputeNodeLoad (dbServiceFactory, loadCase);
//     // }


    
    
//     // auto& assembly = AssemblyEigen::GetInstance();
//     // auto subcaseId = subcase.GetSubcaseId();
//     // assembly.Assembly (dbServiceFactory, subcaseId);
//     // std::cout << "== assembly fini " << std::endl;

//     //std::cout << K << std::endl;


//     return Status::STATUS_Success;
// }


void ThermalSolver::Solve (const Kernel::MModel& model, const Id loadId, _OUT std::vector<SG::DataStructure::Common::Real>& temperature)
{
    auto& assembly = AssemblyEigen::GetInstance();

    EigenVector f;
    model.ComputeLoadVector(loadId, f);

    //std::cout << f << std::endl;
    
    EigenVector temperature0;
    s_solver->solve (f, temperature0);
    if (Eigen::Success != s_solver -> info ())
    {
        std::cerr << "solve failed !!!" << std::endl;
        // return;
    }

    assembly.RecoveryDisplacement (temperature0, _OUT temperature);
    std::cout << "== solve Fini"<< std::endl;
    // 积分点应力应变计算， 外推插值

    // 应力磨平

    std::cout << "Thermal Solve fini  \n";
}


void ThermalSolver::Solve (const DBServiceFactorySP dbServiceFactory, const SG::DataStructure::FEM::Project& project, const SG::DataStructure::Common::IdArray& ids)
{
    // 对应工况的分析类型必须是 `传热`
    // 即 AnalysisTypes::ANALYSISTYPE_Heat == subcase.m_sol
    // AnalysisProcess::Initialize (dbServiceFactory, ids); // 放开函数内部的注释（分析类型）
    heatElementCalculatorRegister ();

    Kernel::MModel model (SOLType::SOL_NLSCSH, dbServiceFactory);
    model.Initialize (ids);
    model.CalcMatrixEquation (ids[0]);
    
    // 初始化稀疏矩阵求解器
    // std::string solver_type = "sgEigenLLT";
    // std::string solver_type = "sgSuperLU";
    // std::string solver_type = "sgPetscKsp";
    std::string solver_type = "TPardiso<Real_t>";//sgPardiso
    auto        solver      = reinterpret_cast<SG::Algebra::LinearSolver*> (SG::Base::Factory::GetInstance ()->Produce (solver_type));
    if (solver)
    {
        s_solver.reset (solver);
    }
    assert (solver && "Make sure that the required solver exists!");
    
    SG::Algebra::SpMat_t K;
    model.ComputeStiffnessMatrix (_OUT K);
    s_solver -> compute(K);

    if (Eigen::Success != s_solver -> info ())
    {
        std::cerr << "Decomposition failed !!!" << "fail type " << s_solver -> info () << std::endl;
        return;
    }
    std::cout << "== Matrix Decomposition" << std::endl;
    
    
    // std::set<Id> loadCases = dbServiceFactory->get<ISubcaseService> ()->GetLoadsetIds (ids);
    // 热分析把热当成载荷
    // for (const auto& subcaseId : ids)
    // {
    //     auto para = dbServiceFactory -> get<ISubcaseService>() -> FindSubCase (subcaseId);
    //     Subcase subcase (para);
    //     auto spcId = subcase.GetSPCId();
    //     if (0 != spcId)
    //     {
    //         loadCases.emplace(spcId);
    //     }
    // }


    // 依次计算每个工况的温度结果
    const auto& temperatureService    = dbServiceFactory->get<ITEMPERATUREService> ();
    const auto& subcaseService  = dbServiceFactory -> get<ISubcaseService>();
    for(auto subcaseId : ids)
    {
        // 获取LOAD的荷载集ID
        auto para = subcaseService -> FindSubCase (subcaseId);
        Subcase subcase(para);
        std::cout << "== Get RHS " << std::endl;
        const auto loadId = subcase.GetLoadId ();   
        // 计算温度
        std::vector<SG::DataStructure::Common::Real> temperature;
        ThermalSolver::Solve (model, loadId, _OUT temperature);

        // 保存subcase工况温度
        std::vector<SG::DataStructure::Post::TEMPERATUREData> temperatureData;
        SG::FEMSolver::DataAdapter::PackGridTemperature (temperature, _OUT temperatureData);
        temperatureService->Add ({SOLType::SOL_NLSCSH, subcaseId}, temperatureData);
    }

    // perf: force to release the resources allocated for the linear solver.
    s_solver.reset ((SG::Algebra::LinearSolver* ) (nullptr));
}